/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of LINAL (LINear ALgebra)
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/
#include "gpu_test.hxx"

int main(int argc, char **argv) {

  linal::Flat_  A1, A2, B, X1, X2, p1, p2, norm;
  linal::Hier_ hA1, hA2, hp1, hp2, hX1, hX2;

  if (argc != 6) {
    printf("Try :: %s ", argv[0]);
    printf("[n_thread] [n_thread_gpu] [n_gpu] ");
    printf("[n_bin] [n_work]");
    return -1;
  }

  int verbose = 0;

  int
    n_thread     = atoi( (argv[1]) ),
    n_thread_gpu = atoi( (argv[2]) ),
    n_gpu        = atoi( (argv[3]) ),
    n_item_bin   = atoi( (argv[4]) ),
    n_item_work  = atoi( (argv[5]) );

  int datatype = TEST_DATATYPE, n_device;

  cudaGetDeviceCount( &n_device );
  // ---------------------------------------
  // ** Initialization 
  FLA_Init();

  // ---------------------------------------
  // ** Matrices
  A1.create   (datatype, N, N);
  A2.create   (datatype, N, N);

  p1.create   (FLA_INT, N, 1);
  p2.create   (FLA_INT, N, 1);

  B.create    (datatype, N, 1);
  X1.create   (datatype, N, 1);
  X2.create   (datatype, N, 1);
  
  norm.create(datatype, 1, 1);

  FLA_Random_spd_matrix( FLA_LOWER_TRIANGULAR, ~A1);
  FLA_Hermitianize( FLA_LOWER_TRIANGULAR, ~A1 );

  FLA_Random_matrix( ~B );

  FLA_Copy(~A1, ~A2);

  hA1.create( A1, GPU_BMN, GPU_BMN );
  hA2.create( A2, GPU_BMN, GPU_BMN );

  hX1.create( X1, GPU_BMN, GPU_BMN );
  hX2.create( X2, GPU_BMN, GPU_BMN );

  hp1.create( p1, GPU_BMN, GPU_BMN );
  hp2.create( p2, GPU_BMN, GPU_BMN );

  if (verbose) {
    FLA_Obj_show("=A1=", ~A1, "%4.3f", "===");
    FLA_Obj_show("=A2=", ~A2, "%4.3f", "===");
  }

  // ---------------------------------------
  // ** Chol incpiv
  omp_set_num_threads( n_thread );

  linal::init_gpu(n_thread_gpu, n_gpu, n_item_bin, n_item_work );
  linal::set_cpu_blocksize(CPU_BMN);
  linal::set_work_balance(1);
  linal::set_computing_model( TEST_COMPUTING );
  linal::create_flat_gpu( datatype, GPU_BMN, GPU_BMN );

#pragma omp parallel
  {
#pragma omp single nowait
    linal::dense::chol_incpiv( FLA_LOWER_TRIANGULAR, hA2, hp2 );
  }
  linal::apply_pivots_ln_internal( ~p2, ~B );

  linal::disp_gpu_data_access(stdout);
  linal::disp_operation_stat(stdout);

  linal::set_computing_model( LINAL_CPU );
  linal::free_flat_gpu();
  linal::finalize_gpu();

  FLA_Chol_solve( FLA_LOWER_TRIANGULAR, ~A2, ~B, ~X2 );

  // ---------------------------------------
  // ** FLAME
  linal::apply_pivots_ln_internal( ~p2, ~A1 );
  linal::apply_pivots_rt_internal( ~p2, ~A1 );

  FLA_Chol( FLA_LOWER_TRIANGULAR, ~A1 );

  FLA_Chol_solve(FLA_LOWER_TRIANGULAR, ~A1, ~B, ~X1);

  if (verbose) {
    FLA_Obj_show("=LINAL X2=", ~X2, "%4.3f", "===");
    FLA_Obj_show("=FLA   X1=", ~X1, "%4.3f", "===");
  }

  // ---------------------------------------
  // ** Check
  FLA_Axpy(FLA_MINUS_ONE, ~X1, ~X2);
  FLA_Norm1( ~X2, ~norm);

  if (verbose)
    FLA_Obj_show("=Diff X2=", ~X2, "%4.3f", "===");

  float  norm_f = *((float *)FLA_Obj_buffer_at_view(~norm));
  double norm_d = *((double*)FLA_Obj_buffer_at_view(~norm));

  int rval;

  printf("- TEST::");
  for (int i=0;i<argc;++i)
    printf(" %s ", argv[i] );
  printf("\n");

  switch (datatype) {
  case LINAL_SINGLE_REAL:
    if (norm_f < LINAL_ERROR_TOL) {
      printf("PASS::Norm :: %E \n", norm_f);   rval = 0;
    } else {
      printf("FAIL::Norm :: %E \n", norm_f);   rval = -1;
    }
    break;
  case LINAL_DOUBLE_REAL:
    if (norm_d < LINAL_ERROR_TOL) {
      printf("PASS::Norm :: %E \n", norm_d);   rval = 0;
    } else {
      printf("FAIL::Norm :: %E \n", norm_d);   rval = -1;
    }
    break;
  }

  // ---------------------------------------
  // ** Matrices
  A1.free(); hA1.free();
  A2.free(); hA2.free();
  p1.free(); hp1.free();
  p2.free(); hp2.free();

  B.free();
  X1.free(); hX1.free();
  X2.free(); hX2.free();

  norm.free();

  // ---------------------------------------
  // ** Finalization   
  FLA_Finalize();

  return rval;
}
